radius=20;
sigma=1.2;
vp=zeros(41,41);
angleImg=zeros(41,41);
for j=1:2*radius+1;
 for k=1:2*radius+1;
    angleImg(j,k)=atan2(k-radius-1,j-radius-1) ;
 end
end

ea=zeros(41,41);
for j=1:41;
 for k=1:41;
%    ea(j,k)=exp(1-(pi/2)/abs(pi/2-abs(angleImg(j,k)-0)+0.0001));
   ea(j,k)=exp(-abs(tan(abs(angleImg(j,k)-0)).^2/(2*0.7*0.7)));
 end
end
% j=1:2*radius+1;
% k=1:2*radius+1;
j=-radius:1:radius;
k=-radius:1:radius;
[xarry,yarry]=meshgrid(j,k);
colormap(gray);
surf(xarry,yarry,ea,ea);



% radius=20;
% sigma=1.2;
% for j=-radius:1:radius;
%  for k=-radius:1:radius;
% [xarry,yarry]=meshgrid(j,k);
% vp=[yarry,xarry];
% sr=xarry.^2+yarry.^2;
% er=exp(-sr/(2*sigma*sigma+0.0001));
% surf(xarry,yarry,er,er);
% ea=exp(1-(pi/2)/abs(pi/2-abs(angleImg(rp(1,1),rp(1,2))-angleImg(vp(1,1),vp(1,2)))+0.0001));



% beta=1;
% j=-radius:1:radius;
% k=-radius:1:radius;
% [xarry,yarry]=meshgrid(j,k); 
% sr=xarry.^2+yarry.^2;
% gx=((Ix(vp(1,1),vp(1,2))-Iy(vp(1,1),vp(1,2))).^2+(2*Ixy(vp(1,1),vp(1,2))).^2)/(Ix(vp(1,1),vp(1,2))+Iy(vp(1,1),vp(1,2))).^2;
% cx=(1-gx)*magImg(vp(1,1),vp(1,2)).^2;
% radiu=sqrt(sr);
% sigma1=radiu/(1+cx/beta);
% sigma2=(1-gx)*sigma1;
% n=[cos(angleImg(vp(1,1),vp(1,2))+pi/2);sin(angleImg(vp(1,1),vp(1,2))+pi/2)];
% nt=[cos(angleImg(vp(1,1),vp(1,2)));sin(angleImg(vp(1,1),vp(1,2)))];
% er=exp(-(((vp-rp)*n).^2/(sigma1^2+0.0001)+((vp-rp)*nt).^2/(sigma2^2+0.0001)));